MODULE md1_hermite_finite_mpoles

   USE md1_global_consts
   USE md1_global_types
   IMPLICIT NONE
   PRIVATE
   ! Public procedures
   PUBLIC :: md1_build_1D_finite_mpoles

CONTAINS

!-------------------------------------------------------------------------------

   SUBROUTINE md1_build_1D_finite_mpoles(xyz,batch,abc,Qjie)
     
      USE md1_hermite_expansion, ONLY: md1_build_Ecoefs
      
      IMPLICIT NONE
      INTEGER(INTK),        INTENT(IN)  :: xyz
      TYPE(md1_prim_batch), INTENT(IN)  :: batch(:)
      REAL(REALK),          INTENT(IN)  :: abc(3)
      REAL(REALK),          INTENT(OUT) :: Qjie(:,0:,0:,0:)

      REAL(REALK), ALLOCATABLE :: Mte(:,:,:)
      REAL(REALK), ALLOCATABLE :: Etij(:,:,:,:)
      INTEGER(INTK) :: np, Iang, Jang, LMAX, IJang

      np    = UBOUND(Qjie,1)
      Jang  = UBOUND(Qjie,2)
      Iang  = UBOUND(Qjie,3)
      LMAX  = UBOUND(Qjie,4)
      IJang = Iang+Jang

      ALLOCATE (Mte(np, -1:Iang+Jang+LMAX+1, -1:LMAX+1))
      ALLOCATE (Etij(np, 0:Iang+Jang, 0:Iang, 0:Jang))
      
      CALL md1_build_hermite_finite_mpoles(batch,abc,xyz,IJang,LMAX,Mte)
      CALL md1_build_Ecoefs(batch,xyz,Iang,Jang,Etij)
      CALL md1_transform_hermite_to_car(LMAX,Etij,Mte,Qjie(:,:,:,:))

      DEALLOCATE (Mte)
      DEALLOCATE (Etij)

   END SUBROUTINE md1_build_1D_finite_mpoles

!-------------------------------------------------------------------------------

   SUBROUTINE md1_build_hermite_finite_mpoles(batch,abc,xyz,IJang,LMAX,Mte)

      ! Evaluate 1D hermite multipole integrals on the interval [a,b] about c

      IMPLICIT NONE
      TYPE(md1_prim_batch), INTENT(IN)  :: batch(:)
      REAL(REALK),          INTENT(IN)  :: abc(3)
      INTEGER(INTK),        INTENT(IN)  :: xyz, IJang, LMAX
      REAL(REALK),          INTENT(OUT) :: Mte(:,-1:,-1:)

      REAL(REALK) :: TA(SIZE(batch), 0:LMAX+IJang)
      REAL(REALK) :: TB(SIZE(batch), 0:LMAX+IJang)
      REAL(REALK) :: expI(SIZE(batch))
      REAL(REALK) :: PC(SIZE(batch))
      REAL(REALK) :: a(SIZE(batch))
      REAL(REALK) :: b(SIZE(batch)), tt
      INTEGER(INTK) :: e, t

      ! Set integration range [a,b]
      a(:) = batch(:)%P(xyz) - abc(1)
      b(:) = batch(:)%P(xyz) - abc(2)
      ! Set multipole origin c
      PC(:) = batch(:)%P(xyz) - abc(3)

      Mte(:,:,:) = zero
      TA(:,:) = zero
      TB(:,:) = zero

      ! t=0; e=0 
      TA(:,0) = ERF( SQRT(batch(:)%exp) * a(:) )
      TB(:,0) = ERF( SQRT(batch(:)%exp) * b(:) )
      Mte(:,0,0) = half * SQRT(PI/batch(:)%exp) * (TA(:,0) - TB(:,0))

      ! t=1; e=0 
      TA(:,1) = EXP( -batch(:)%exp * a(:)*a(:) )
      TB(:,1) = EXP( -batch(:)%exp * b(:)*b(:) )
      Mte(:,1,0) = TA(:,1) - TB(:,1)

      ! t=2..(LMAX + IJang) ; e=0 
      DO t = 2, LMAX + IJang
         tt = DBLE(t)
         TA(:,t) = -two * batch(:)%exp * (a(:)*TA(:,t-1) + (tt-two)*TA(:,t-2))
         TB(:,t) = -two * batch(:)%exp * (b(:)*TB(:,t-1) + (tt-two)*TB(:,t-2))
         Mte(:,t,0) = TA(:,t) - TB(:,t)
      END DO

      ! Recurrences for remaining {t,e} same as usual, with different t-range
      expI(:) = half / batch(:)%exp
      DO e = 1, LMAX
         DO t = 0, LMAX + IJang - e
            Mte(:,t,e) =  PC(:) * Mte(:,t,e-1)           &
                      + expI(:) * Mte(:,t+1,e-1)         &
                      + DBLE(t) * Mte(:,t-1,e-1)
         END DO
      END DO

   END SUBROUTINE md1_build_hermite_finite_mpoles

!-------------------------------------------------------------------------------

   SUBROUTINE md1_transform_hermite_to_car(LMAX,Etij,Mte,Qjie)

      IMPLICIT NONE
      INTEGER(INTK), INTENT(IN)  :: LMAX
      REAL(REALK),   INTENT(IN)  :: Etij(:,0:,0:,0:)
      REAL(REALK),   INTENT(IN)  :: Mte(:,-1:,-1:)
      REAL(REALK),   INTENT(OUT) :: Qjie(:,0:,0:,0:)

      INTEGER(INTK) :: i,j, e,t

      Qjie(:,:,:,:) = zero

      DO e = 0, LMAX
         DO i = 0, UBOUND(Etij,3)
            DO j = 0, UBOUND(Etij,4)
               DO t = 0, i+j
                  Qjie(:,j,i,e) = Qjie(:,j,i,e) + Etij(:,t,i,j) * Mte(:,t,e)
               END DO
            END DO
         END DO
      END DO

   END SUBROUTINE md1_transform_hermite_to_car

!-------------------------------------------------------------------------------

END MODULE md1_hermite_finite_mpoles
